function [CI_logQ_low, CI_logQ_high]=compute_CIs_byincome(qhat,qhat_star,reps,p_interest,y_interest,p_interest_univ,y_interest_univ)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % INITIALIZE
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

points_allincomes=size(p_interest,1);
points_byincome=size(p_interest_univ,1);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE JOINT STANDARD ERRORS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CI_logQ_low=NaN*ones(3,points_allincomes);
CI_logQ_high=CI_logQ_low;
alpha=NaN*ones(size(y_interest_univ,1),1);

for yy=1:size(y_interest_univ,1),

    select=(y_interest==y_interest_univ(yy,1));
    fprintf('now joint CI for yy=%6.3f\n',yy);
    [CI_logQ_low(2,select), CI_logQ_high(2,select),alpha(yy,1)] = find_jointCI(qhat(select,1),qhat_star(select,:)',points_byincome,p_interest(select)',reps);

end;

    %
